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I. INTRODUCTION 



00 

The cosmological singularity in spatially homogeneous models is known to be either asymptotically velocity term 
dominated (AVTD) (Kasner-like) gj|| or Mixmaster-like §g]. Analytic |J| and numerical studies (|J^] have shown 
, that the singularity in the spatially inhomogeneous Gowdy models || is AVTD everywhere in polarized models and 
■ everywhere except, perhaps, at a set of measure zero in the generic, unpolarized case. Long ago, Bclinskii, Khalatnikov, 
and Lifshitz (BKL) || argued that the generic singularity in Einstein's equations is locally of the Mixmaster type. 
t-H i This remains controversial Q but serves as a conjecture which can be tested. The Gowdy models, though interesting 
in their own right, have two commuting Killing fields which precludes local Mixmaster behavior for T 3 x R topology. 
(For generalized Gowdy models which appear to exhibit Mixmaster dynamics see [p~0| .) However, the more general 
one Killing field t/(l) symmetric cosmologies on T 3 x R should gcncrically have Mixmaster- like singularities if the 
BKL conjecture is correct. 

An AVTD solution to Einstein's equations is obtained by neglecting all terms containing spatial derivatives. The 
parameters of the solution to the resultant ODE's are then assumed to depend on the spatial coordinates. An AVTD 
0^ ■ singularity is then one where the solution to the full Einstein equations comes arbitrarily close to an AVTD solution as 
*q \ the singularity is approached For the U(\) cosmologies, the AVTD equations may be solved exactly. However, one 
QT* may worry about the coordinate and slicing dependence of the notions of space and time leading to such dependence 
in our definition of AVTD. We shall argue here that the AVTD behavior of the singularity is not strongly slicing 
(3J[). dependent. 

We have already shown || that numerical methods used in the Gowdy cosmologies can be easily generalized to the 
. j—{ \ U(l) case. However, numerical difficulties associated with spatial differencing in two dimensions have prevented so 
far a complete understanding of generic U(l) models probably because Gowdy- like spiky features Q cannot yet be 
modeled with sufficient accuracy. However, the subclass of polarized C7(l) models does not appear to develop spiky 
features and thus can be treated well by our numerical methods. We consider vacuum £7(1) symmetric cosmologies 
on T 3 x R. These can be described by two variables ip and u> analogous to the Gowdy wave amplitudes for the + 
and x polarizations of gravitational waves and three variables A, x, and z which contain nondynamical information. 
Polarized U (1) models have the degree of freedom associated with ui set equal to zero. This condition is preserved for all 
time both by Einstein's equations and by the discrete form thereof used in our numerical simulations. While a general 
solution to the initial value problem is not known explicitly, we have found an algebraic solution to the constraints 
which contains three (four for generic U(l) models) arbitrary functions of the spatial variables and is sufficiently 
general to yield generic evolution toward the singularity within the polarized U(l) class of models. We find that the 
evolution toward the singularity is AVTD everywhere as had previously been conjectured This conclusion is 

based (a) on the exponential decay with time of all terms in the Hamiltonian containing spatial derivatives, (b) the 
exponential decay with time of the change in certain functions of the spatial coordinates which are constant in time 
in the AVTD regime, and (c) fits of the behavior of the variables at randomly selected spatial points to the expected 
linear or constant asymptotic time dependence. 
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II. THE MODEL 



The spacelike Killing vector £ = d/dx 3 of the U(l) symmetry allows a generic metric in this class to be described 
in terms of the coordinate x 3 along £ and coordinates u and v and r parametrizing the quotient 2 + 1 manifold [fL^ . 
Using notation from f| and taking u 6 [0, 27r] and v <E [0, 2tt] the 3 + 1 metric is 

ds 2 = e~ 2ip {-N 2 e~ 4T dT 2 + e- 2T e A e ab dx a dx b } + e 2v (dx 3 + a dx a ) 2 (1) 
where {x a } — {u 7 v} = {x 1 ,^ 2 }, the lapse N = e A , and g ab = e A e a b where 



e ab - 2 



e 2z + e- 2z (l + x) 2 e 2z + e- 2z (x 2 - 1) 
e 2z + e- 2z (x 2 - 1) e 2z + e~ 22 (l -x) 2 



(2) 



is the conformal metric of the u-v plane. 

In the polarized case, (3 a = but for the general case we must construct {/3 a } = {Pitfo} as follows: At an initial 
time, t = t , define 

(3$(u,v) = ci + / dv'\r{u,v')-^- [ dv"[r(u,v")] 
Jo I 2tt J 

pUu,v) = c 2 - — / du' / dv'r(u>') (3) 
27r Jo Jo 

where ci and C2 are constants. These constants are like "twist" constants [|l2| and are arbitrary (but physically 
significant). These (3*'s will be periodic on T 2 provided 

2tt p2-rr 

du / dv r(u,v) = 0. (4) 
Jo 

Now let pi = fli + dX/du and = (3% + d\/dv where A is arbitrary and could be identically zero to give {(3 a } at t . 
Finally, define (3 a (u,v,T) by computing 

f3 a (u, v, t) = a (u, v, to) - f dr' e - 2T ' Ne-^e ab s bc u, c (5) 

where 

_ 6c ( o :i 



£ -^-1 oj- ( 6 ) 

As in [p^| , a canonical transformation replaces the terms in the Einstein-Hilbert action e a (3 a + 0o e a , a by rui where 
the momentum e a conjugate to j3 a identically satisfies the constraint e a , a = if e a = e a u>, b . Einstein's equations 
may be obtained by the variation of 



H = J J dudvH 

= JJ dudv Qp 2 + \^ z V l + \ V 2 + \e^r 2 - \p\ + 2p A ) 

+e- 2T f fdudv{ (eV b ) , ab - (eV 6 ) , a A, b +e A [(e~ 2z ) ,„ x, v - {e~ 2z ) , v x, u ' 



+2e W,a <p, b +^e A e-^e ab io, a uj, b 

= H K + H V = J JdudvH K + J J dudv V. (7) 

The Hamiltonian and momentum constraints are respectively 

H° = U - 2 PA = (8) 
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and 



H u = p z z, u +p x x, u +pa A,„ -pa,u +Pf,u +rcj, u +- { [e iz ~ (1 + x) 2 ] p x - (1 + x)p z ) , v 
-\{[e iz + {l-x 2 )] p s -xp z }, u =0, 



(9) 



' Pz Z-iv \~Px 



"PA, 



4z 



{\~x) 2 ]p x + {\-x)p z } , u 



+ 2 {l e4Z + i 1 ~ x 2 )} Px~xp z }, v = 0. 



(10) 



To evolve Einstein's equations numerically, we need to solve the initial value problem. While an explicit general 
solution is not available, a particular class of solutions may be obtained as follows: To solve the momentum constraints 
(0) and ( |To| ) set p x — p z — ip, a — u, a — to leave pa A, a —pA,a = which may be satisfied by requiring pa — ce A . For 
sufficiently large c, the Hamiltonian constraint may be solved algebraically for either p or r. In general, this leaves as 
free data the four functions x, z, A, and either r or p. However, in the polarized case, we demand 



u) = r = 



(11) 



so that the Hamiltonian constraint must be solved for p with three arbitrary functions x 7 z, and A. We note that 
while any c will solve the initial value problem, to approach the singularity we require c > in order to give pa > 
so that A will decrease (become more negative) as t — > oo. This is required because the determinant of the metric 
has a factor e 2A which measures the area in the u-v plane. 

The AVTD equations may be obtained by variation of Hk in (0) and have the exact solution (see PJl^|) 



z = -v z (t - t 0z ) + ln[\fi z \ (1 + e - 4v ^-^)] - Vl 



t as t 



oo, 



X 




Mi + 


e -iv z {T-T z)\ 


-1 


xq as r — > oo, 


Pz 


= -4w 2 


(l-e- 






— Av z as T — > OO , 


(1 + e- 


4:Vz(t-T z)\ 




Px 


= -n z v z =p° x , 










= -«<p( T - T 0v 


+ HM0- + 




)] — > — w^r as r — > oo 


U) 




H v (1 - 




-1 


-> wo as t — > oo, 


P 




(1-e- 






-> — 4u y as r — > oo, 


(1 + e- 







r = -p, v v v = r , 
A = A + (2-p°)r, 

PA =P°A 

subject to the AVTD limit of the Hamiltonian constraint (g) 



pi 



i 



i 



-pi + e* z p£ + -jf + e^r 



<t<P r 2 



(12) 



(13) 



The AVTD limit is expressed in terms of \i z , v z > 0, £ z , T 0z, ^ip, v v > 0, and tq v which are functions of u and v. 
The limits on v z and v v arise in order for the limiting forms of (12) to cause the ex pon ents in (0) to decay as r — > oo. 
(See the discussion of similar terms in Gowdy models HQ].) The AVTD solution ([12]) may be inverted to give these 
parameters in terms of the original variables. Useful examples are 



Y^P 2 Z + \^ z Pl 



—p 2 + -e^r 2 
16 P 4 



We further note [^J that, in addition to v z , v v , p x , r, and pa, 

1 



2 Pz 



PxX 



(14) 



(15) 
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are also constant in time in an AVTD regime. 

The method of Grubisic and Moncrief (GM) § can be followed to determine if one expects the singularity to 
be AVTD in the polarized model. We req uire every term in V (from (Q)) to decay exponentially as r — > oo if we 
substitute the limiting AVTD solution, (|l2|), for the U(l) variables. Clearly, the behavior will depend on the behavior 
of the exponentials. First, as in the Gowdy case we expect the remaining nonlinear term in Hk, s 4z p 2 ,/2 to 

drive p z to negative values yielding z eventually large and negative. In that case, terms in e ab proportional to e~ 2z 
will dominate. From (]?]), we must examine the factor e -2r+A-2z wn i cn w jH decay as r — > oo if p\ — 2v z > 0. (Recall 
that v z > 0.) But the AVTD limiting form of the Hamiltonian constraint requires 

p 2 + p 2 - 4p| = = 16v 2 +p 2 - Ap\ (16) 

which, in turn, requires pa > 2v z . Thus, p\ > is expected if the singularity is AVTD and yields consistent behavior. 



III. NUMERICAL METHODS 

To integrate Einstein's equations, the symplectic method described in detail in jf| is applied to U(l) symmetric 
cosmologies. For a system described by a Hamiltonian 

H = H 1 + H 2 , (17) 

the corresponding evolution operator U(At) which evolves data from r to r + At generated by H can be written as 

U(Ar) = [ / 1 (^)f/ 2 (Ar)C7 1 (^) + G[(Ar) 3 ] (18) 

where U\,2 are the evolution operators generated by i?i,2- Suzuki |l6|Jrif | has shown how to generalize this algorithm to 
an arbitrary order in Ar. This method is useful if the sets of equations of motion obtained by the separate variations 
of Hi and if 2 are exactly solvable. 

We see Q that ([?]) for the U(l) models has exactly solvable equations for Hk and Hy- We have already obtained 
the AVTD solution ( |l2] ) from Hk's equations of motion. Since Hy contains no momenta, the configuration variables 
x, z, A, tp, and u are constants of the motion so the equations are trivially solved. Non-trivial, however, is the 
representation of the gradients of V which arise during the variation. Accurate representation of spatial derivatives in 
two spatial dimensions is more difficult than in one dimension since the function of interest is no longer guaranteed to 
vary along the differencing direction. These inaccuracies have been found to limit the duration of generic U(l) model 
simulations but are not problematical for polarized models. This is almost certainly due to the fact that Gowdy-like, 
steepening spiky features characterized by increasingly (with r) large spatial gradients are absent in polarized U(l) 
models as they also are in polarized Gowdy models. We emphasize that the "averaging" mentioned in [1^] as necessary 
to stabilize generic U(l) evolution is not required in the polarized case. 

To illustrate the behavior of typical polarized U(l) symmetric models, we consider a particular choice of initial 
data with A = Asinusinv, x = z = cosmcosi>, and p\ = ce A with A = 1, c = 14. Other choices of initial data yield 
qualitatively similar results. 

The constraints are imposed initially and monitored thereafter. Fig. 1 shows the maximum value of the constraints 
on the spatial grid vs time and spatial resolution. Although it is not clear how small one should require the constraints 
to be, the demonstrated convergence allows one to argue that the small constraint violation can be neglected in our 
interpretation of the results. Choptuik has argued that constraint convergence is an indicator of convergence of the 
numerical solution to the true solution to Einstein's equations [fL8[ . 



IV. RESULTS 

Fig. 2 displays the maximum value of log 10 V from (^) at two spatial points as a function of r. One expects this 
quantity to vanish as r — > 00 for an AVTD singularity. Here we demonstrate that it vanishes exponentially as one 
expects. We note that the exponential decay of V measures only consistency and does not prove AVTD behavior. 
The observed exponential decay is compared to the predicted e^ _PA+Pz / 2 ' T behavior. While only representative spatial 
points have been shown, exponential decay at the predicted rate is seen at all spatial points. 

In Fig. 3, the maximum values over the spatial grid of the change with time in the AVTD functions that should be 
constant in time in the AVTD limit — p x , v z , v^, c z , c v , pa — are plotted vs r. Again, one sees exponential decay to 
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the level of machine noise. (This machine noise does not show up in log 10 V which is an evaluation but does appear 
when differences are computed.) These quantities are strictly constant in the AVTD limit. In order to compare the 
observed to predicted decay of these quantities, one needs to go beyond the AVTD solution. In the Gowdy case, this 
was done by GM § . 

Finally, in Fig. 4, the variables A, z, and <f>, expected to be linear with r in the AVTD limit, and x expected to be 
constant in r are shown vs r at typical spatial points. For the first three, linear fits are shown. We see that x does 
not show the constant in r behavior one would expect as r — > oo. However, we note from Fig. 4 that x w 0. This 
small constant value allows the exponentially decaying term in x to be measurable. The lines in Fig. 4(d) are fit by a 
constant plus an exponential decaying with the expected rate of Av z from (|l2|). Again, although only representative 
points are displayed, the same behavior is seen at all spatial points. 

V. DISCUSSION 

Numerical studies of polarized U{\) models provide strong evidence in support of their conjectured AVTD singular- 
ity. For a representative choice of initial data, all relevant terms are examined for consistency with AVTD behavior. 
Terms predicted to decay exponentially in the AVTD limit as r — > oo do so with the expected slope. The terms c z , 
Cy, v z , Vip, p x , and p\ are strictly constant in r for all r in the AVTD solution. Decay to the observed constant values 
could be compared with the prediction from the next order in an asymptotic expansion where the AVTD limit is the 
zero order term [^[. This has not been done. However, this is not required to demonstrate that the singularity is 
AVTD. 

Given the consistency of our numerical results with the AVTD limit, we argue that longer duration simulations and 
other initial data sets will reveal nothing new. This is because the method of GM suggests that exponential decay of 
all the terms in V is consistent with the AVTD limit as r — > oo. As shown in the Gowdy case || and, of course, in 
Mixmaster itself || , the exponential growth of such terms is required to change the qualitative behavior of the model 
at a given spatial point. In the absence of this growth, there is no mechanism within Einstein's equations that could 
destroy the AVTD behavior. 

Finally, we consider the possible slicing dependence of the notion of AVTD. First, we note that the canonical 
transformation to r and to has interchanged generalized coordinate and conjugate momentum. Of course, this degree 
of freedom is absent in polarized U(l) models. However the algebraic coordinate conditions which we have imposed in 
fixing the metric form (1) (essentially zero shift and a "harmonic" time function when the choice N = e A is enforced) 
are not completely rigid since they depend still upon the choice of an initial hypersurface. There are many other 
harmonic time functions (i.e., solutions of the wave equation having timelike gradient) which could have been used 
instead of any given one without disturbing the metric form but one can analyze these asymptotically using the formal 
expansion methods of GM. This (not yet rigorously justified) analysis determines the asymptotic behavior of any such 
time function relative to a given one (with respect to which AVTD behavior of the metric is assumed) and shows that 
asymptotic velocity term dominance would be preserved in any harmonic time, zero shift gauge provided it holds in 
the original one. Thus the occurence of AVTD behavior does not seem to depend upon the choice of a preferred initial 
hypersurface, at least not within the class of gauges under study. 
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FIGURE CAPTIONS 



Figure 1. Convergence of the constraints. The Hamiltonian (triangles), u-momentum (circles), and v- momentum 
constraints vs r for spatial resolutions of 128 x 128 (broken line), 248 x 248 (solid line), and 512 x 512 (dashed line), 
(a) A second order (in time) accurate symplectic PDE solver yields a decrease in the magnitude of the constraints 
by a factor of 4 as the number of spatial grid points in each direction is doubled (second order convergence), (b) 
A fourth order (in time) accurate symplectic PDE solver jl6],[l7) yields a factor of 16 decrease in the magnitude of 
the constraints (fourth order convergence). Only the early part of the evolution is shown. Note that in addition to 
improvement in convergence with increasing spatial resolution, one also finds convergence with increasing time order 
accuracy. In all cases, the spatial differencing scheme |l9) is fourth order accurate and centered about the spatial grid 
point of interest . 

Figure 2. Exponential decay of the spatial derivative terms in the Hamiltonian. log 10 \ V\ vs r is displayed for two 
spatial points (solid line). The straight line slope indicates exponential decay. In each case, the corresponding values 
of (— pa +Pz/2)r/lnlO are also displayed (dashed line). Figures 2-5 are obtained from a second order accurate PDE 
solver with 512 2 spatial grid points. 

Figure 3. Decay of the deviations from constant in time behavior of quantities constant in the AVTD limit. The 
maximum value of the change with time over the spatial grid of the quantities p x , c z , v z , c v , v v , and p\ are shown vs 
r. The non-vanishing values at late times are at the level of machine precision errors. 

Figure 4. Comparison of predicted and measured values of the variables (a) A, (b) z, (c) tp, and (d) x vs r are 
shown for two representative spatial points. For (a), (b), and (c), the data is represented by circles or squares while 



the solid line is a linear fit to the data. For (d), the fit is to x — a + be P * T as is consistent with (12) 
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